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ABSTRACT 

Numerical solution to a theoretical model of vapor cavitation in a 
dynamically loaded journal bearing Is developed, utilizing a multigrid 
iterative technique. The method Is compared with a noniterative approach in 
terms of computational time and accuracy. The computational model is based on 
the Elrod algorithm, a control volume approach to the Reynolds equation which 
mimics the Jakobsson-Floberg and Olsson cavitation theory. Besides accounting 
for a moving cavitation boundary and conservation of mass at the boundary, it 
also conserves mass within the cavitated region via a smeared mass or striated 
flow extending to both surfaces in the film gap. The mixed nature of the 
equations (parabolic in the full film zone and hyperbolic in the cavitated 
zone) coupled with the dynamic aspects of the problem create interesting 
difficulties for the present solution approach. Emphasis is placed on the 
methods found to eliminate solution instabilities. Excellent results are 
obtained for both accuracy and reduction of computational time. 

NOMENCLATURE 

AR aspect ratio of grid size, Ax/A z 
ed dynamic eccentricity, m 



e s static eccentricity, m 

F k forcing function on grid k 

fk residual function on grid k 

g switching function 

H dimensionless film thickness, h/AR 

h film thickness , m 

k-1 

1^ interpolation function from grid k to grid k-1 

k grid indicator, k < M 

L {©} differencing scheme acting on a variable 0 
L/D length-to-dlameter ratio 

M represents the finest grid (highest number) 

MJ number of grid points axially 

m x ,m z lineal mass flux, kg/m-s 
NI number of grid points circumferentially 

p fluid pressure, N/m 2 

p a ambient pressure, N/m 2 

p c cavitation pressure, N/m 2 

R radius of journal , m 

t time, s 

U sum of the surface velocities in x-directlon, m/s 

V sum of the surface velocity vectors, m/s 

WU work units 

x coordinate along circumference, m 

y coordinate normal to x,z plane, m 

z axial coordinate, m 

(3 liquid bulk modulus, N/m 2 

Y angular position of minimum film, rad 
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AR radial clearance, m 

At time Increment, s 

Ax Incremental spacing along circumference, m 

Az axial Incremental spacing, m 

c eccentricity ratio 

0 fractional film content In cavltated zone 

density ratio, p/p c , In full film zone 

\i dynamic viscosity, N-s/m 2 

p fluid densl ty, kg/m^ 

p c fluid density within cavltated zone, kg/m^ 

<p angular coordinate along circumference, degree 

ujj orbital angular velocity of journal center about a fixed point 
relative to the housing center, rad/s 

u>s angular velocity of journal about Its own center, rad/s 

INTRODUCTION 

The presence of vapor cavitation In dynamically loaded journal bearings 
has become a topic of Increasing Importance. The use of Increased loads and 
more complicated loading cycles has resulted In an increase In the occurrence 
of cavitation erosion problems. Examples of journal bearing applications 
Include main and crankshaft bearings In diesel engines and a variety of 
bearings In the aircraft Industry, Dowson and Taylor U). Dynamic loading can 
also lead to Instabilities In the motion, such as whirling or whipping motion, 
which may damage the bearing. In order to avoid bearing damage, It is useful 
to predict the conditions under which the bearing will remain stable. The 
determination of these stability maps requires a knowledge of the hydrodynamic 
force terms. 
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Previous static loading models, such as those using Swift-Stleber or 
Giimbel boundary conditions, assume a stationary cavitation bubble and are 
inadequate for high speed, dynamic applications. Under dynamic loading, 
changes in the local film thickness cause the bubble to grow, move downstream 
from the minimum film position, and collapse, Brewe (2). 

A film model which effectively deals with dynamic loading has been 
formulated by Jakobsson-Floberg (3) and Olsson (4). Besides accounting for a 
moving boundary, it also accommodates the flow within the cavitated region, 
manifested via a smeared mass or liquid striations. These striations have been 
observed In past experimental work. Because the theory assumes a zero pressure 
gradient within the cavitated zone, this mass flow is a Couette flow. The JFO 
theory accounts for both film rupture and film reformation, another advantage 
over previous methods. Unfortunately, the complexity of the JFO theory makes 
it difficult to apply (2). Elrod and Adams (5,6) have developed an algorithm 
which automatically conforms to the JFO theory while being much simpler to 
code. It utilizes a switching function which eliminates the pressure gradient 
terms from the Reynold’s lubrication equation at cavitated points. 

A solution to the Elrod algorithm for a dynamically loaded problem has 
been formulated by Brewe (2). This direct solution to the finite differenced 
equations, i.e., no iterations required, utilizes an alternating direction 
implicit (ADI) scheme for the time march. When compared to a nonconservative 
film model (pseudo-Gumbel boundary conditions), as much as a 20-percent 
difference in load capacity is observed. Brewe's results agree excellently 
with the experimental work of Jakobsson and Floberg (3) for stationary 
cavitation. The present experimental data on nonstationary cavitation is 
limited, but Brewe's results compare reasonably well with the experimental 
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work of Jacobson and Hamrock, (7,8). Unfortunately, Brewe's direct method 
requires two to three times the computational work needed by the nonconservative 
solution. The practical use of the Elrod algorithm In solving dynamic loading 
problems In Industry requires the development of a more efficient computer 
solution. 

Iterative techniques are not considered among the fastest methods of 
computer solution. However, newly developed techniques using multiple grid 
sizes have shown that It Is possible to greatly reduce computational time. It 
Is the purpose of this work to Implement the multigrid technique developed by 
Ach 1 Brandt, (9,10), In the Iterative solution of the Elrod algorithm under 
dynamic loading conditions. As a point of reference, the multigrid method has 
also recently been Implemented In an EHD lubrication problem by Lubrecht, 
ten Napel, and Bosma, (]_]_). 

The nature of the problem, l.e., the presence of discontinuous 
coefficients In the cavltated region, poses Interesting difficulties in the 
application of the multigrid method. Therefore, one major objective Is to find 
methods of making the multigrid technique as effective over the cavltated area 
as It Is over areas of full film. 

DESCRIPTION OF THE MODEL PROBLEM 

The bearing motion consists of a journal undergoing a constant rpm as well 
as a noncentered circular whirl Inside of a 360° cylindrical bearing. The 
journal center moves through a prescribed dynamic cycle, (Fig. 1) from a 
minimum through a maximum eccentricity (see Table I for operating conditions). 

The theoretical model assumes conditions of heavy loading, l.e., load 
carrying capacity >> surface tension forces In the liquid. An oil lubricant 
is used for which the vapor pressure is very nearly zero. It is assumed that 
the lubricant has been degassed and that only vapor cavitation is present. In 
actual experimental work with submerged bearings, this is a necessary 
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condition, since gaseous cavitation forms at near ambient pressure and might 
prevent the occurrence of subambient pressures. The flow balance Into and out 
of the bearing cannot be maintained In the absence of subambient pressures. 
Within the cavltated region, a zero pressure gradient Is assumed. In order to 
determine the load carrying capacity of the bearing under the prescribed 
motion, the flow equations are solved for the pressure variable. It is more 
convenient, however, to introduce another dependent variable, 6, which has a 
dual Interpretation regarding the full film region and the cavitated region. 

In the full film (© > 1.0), © Is the ratio p/p c , which represents the ratio 
of film mass content to the mass content that would exist at the cavitation 
pressure p c . In the cavitated region (© < 1.0), 9 represents the mass 
content that exists in the form of liquid striatlons. Pressure and density 
are related through the bulk modulus, 3, such that (p3p/3p) = 3. A switching 
function (g(©)) automatically eliminates the pressure term at cavitated points 
of the flow. That Is, 

Uncavitated point (© > 1.0): g = 1 

Cavitated point (© < 1.0): g = 0 

THE GOVERNING EQUATION AND DIFFERENCING SCHEME 

The Reynold's lubrication equation written In terms of the fractional film 
content, the switch function, and the bulk modulus becomes 

+ (I) * ^ (eh) - $ • ^j^g ( ©)^ e 

Note that the right-hand side, the ‘pressure Induced flow, completely 
disappears in the cavitated region where the switch function becomes zero. The 
finite difference equations are obtained using a control volume approach as was 
used by Brewe (2). A control volume (Fig. 2) Is constructed about each nodal 
point and the net change In mass flow into the cell (m x+ ax/2 ~ | i 1 x-Ax/2 ) Is 
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equated to the total Increase of mass (p c 0h> in the cell over a time At. The 
control volume equivalent to the mass conservation equation is 



In the full film region (all g = 1), both the convective and pressure 
terms are central differenced, appropriate for the parabolic system. In the 
fully cavitated region (all g = 0) , the pressure terms are eliminated and the 
convective terms are upwind differenced, 9 accounting for the mass transport. 
The combination of switching terms at the cavitation boundary automatically 
sets well posed boundary conditions between the two systems. The resulting Am 
terms from the control volume analysis are 

(“x)c 0 n V ' p c(f) j h -lO * 9 -lXl * h o0 " 9 o)®0 * 1 [ 9 -l h -l( 2 ' 9 o) 

* 9 o h o( 9 -i - 2 * 9 l) - 9 i h i 9 o]} 

Wpress ' (S)(h) f h -l/2 9 -l( e -l - 0 - 0-1/2 * h ?/2>o0o ' 0 

- [ h ?/2 9 i(°i - 03j 

An Euler implicit time differencing scheme is used for stability purposes, 
giving: 

9h - e*h* A % Aril z 
At " p c Ax + P(; Az 

where 0*h* signifies time t - 1. 

It should be noted that all terms on the right-hand side of this equation 
are evaluated at time t and are therefore unknown. Along the axial 
boundaries, l.e., along the edge of the bearing, the boundary condition is 
that of atmospheric pressure. The circumferential direction has wrap around 
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boundary conditions. The problem and boundary conditions are described In more 
detail In (2). 

THE MULTIGRID METHOD 

Analysis of the Single Grid and the Residual Function 

The Iterative solution of a set of equations on a single grid generally 
has rapid convergence over the first few sweeps, but very slow convergence over 
most of the process. By examining the solution process, the reasons for this 
become clear. Assume the following continuous differential equation, 

L{©( x ) > = F(x); with suitable boundary conditions 
for which the discretized set of equations on one grid take the form 

L k{ e k} „ pk (1) 

In the present notation, k represents the particular grid size used, 
represents the differencing scheme acting on 0, and 0 is the exact solution 
of the differenced equations. 

Let 0 be the present approximation to the exact solution 0. By 

substituting 0 Into the differenced equations, the following Is obtained. 

f k = F k _ L k{Qk} (2 ) 

where fk is referred to as the "residual function." The residual function 
Is a means of analyzing the error left in the present approximation. 

A Fourier analysis of fk breaks the error into Its high- and low- 
frequency components. High frequency is defined as wavelength less than or 
equal to four times the grid spacing. After a few relaxation sweeps, e.g., 
Gauss-Seidel , these high-frequency components are smoothed out, due to the fact 
that they are locally corrected. Once the error is all low frequency, the 
smoothing rate drops drastically. The grid spacing is too fine to efficiently 
smooth these low-frequency terms Brandt (9). 
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The Roles of High- and Low-Frequency Error Smoothln 

The basic thrust of multigrid Is to utilize coarser grids to handle these 
low-frequency error terms. As long as the residual function, i.e., the error. 
Is well represented on a particular grid, that grid can quickly smooth its own 
"high-frequency" terms and send an appropriate correction back to the finer 
grid. By utilizing coarser and coarser grids, all of the low frequency terms 
are treated similarly. Besides their ability to efficiently deal with low- 
frequency terms, coarse grids also have fewer nodal points to sweep through, 
making the coarse grid sweeps very cheap. 

Whereas moving to coarser grids smooths the low-frequency error, fine grid 
updates during the multigrid process have value in Improving the accuracy of 
the present solution. The role of relaxation on the finer grid is to resolve 
its own high-frequency components as well as smooth the high-frequency error 
which Is produced by Interpolation from the coarser grid. 


Coarse Grid Representation 


The fine grid problem Itself, I.e., L^t©^} = F k , Is not what Is really 


being represented on the coarse grid k-1 . The actual purpose of the coarse 

k-l 

grid Is to solve for a ’correction value, e£ 


e k_1 = e k ' 1 - e k_1 
c 


as a function of the amount of residual error that is left in the present 
approximation on the fine grid, I.e., f k . 


For a nonlinear problem, the full approximation storage (FAS) mode must be 

used. This method stores the entire value of 9 on the coarse grid k-l 

k-l 

instead of just the correction value © c . If used on a linear problem, the 
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equations reduce to those for the linear mode Brandt (9). The existence of 

nonllnearltles In the convective mass flow contribution necessitates the use 

of the FAS mode. Using FAS, the coarse grid problem becomes 

S *-1 A-l/W 

initial ‘ *k V® ) 

k-1 

where 1^ Is the interpolation operator from the fine to the coarse grid. 
A1 so, 

F k_1 = I k_ 1 (F k - L k e k ) + L k_1 {e k-1 } 

Solve 


L k_1 {e k 1 } = F k_1 


(3) 


ft k-l gk-1 5^-1 
e c " 0 " e In1t1al 


Using the interpolation operator 1 . e . , from the coarse to the fine 

grid, we obtain 

rk 


®new ■ § old * 


k-1 


(4) 

k 


The fact that the coarse grid is solving for 0 C ~ as a function of f 

has two Important consequences. The first Is that f must be well represented 

k-1 

on the coarse grid In order that © c is an accurate correction value. 

k 

Therefore, care must be taken that f Is well smoothed on the fine grid before 
transferring It to the coarse grid. 

The second consequence Is that. It is not necessary that the coarse grid 

u i -k-1 

differencing scheme, L {0 }, exactly match the fine grid differencing 

k ~k 

scheme, L {0 }. This can be seen from the following consideration, I.e., the 
concept of solving for a correction value on the coarse grid arises because 


k-1 r^k-1 


{0"-'} - L 




should 

represent 


[l K { 0 k } - L k {0 k }]. 
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It can be shown that the multigrid equations follow from this basic 
concept. Now considering the RHS, 6 k Is unknown, although we can use the 
fact that l>{e k } E F k . Therefore, 

RHS = [F k - L k {0 k }] = f k = Residual function. 

~k-l k-1 

From the definition of 0 , and the fact that L Is a linear operator, 

V 

the LHS can be written 

L k_1 {® k_1 } = L k “ 1 {6 k " 1 } - L k “^ {0 k_ ^ } 

The full equation then becomes the standard linear multigrid equation 

which can then be adapted to the nonlinear form expressed In Eq. (3). 

Therefore, the terms which must closely represent each other are the 
bracketed terms, not the Individual components within. Thus there Is some 
flexibility In creating L k-1 {e k_1 }. Especially if the problem contains 
rapidly changing spatial coefficients, the coarse grid differencing scheme will 
have to be a modified version of the fine grid scheme. 

Multlqrld Cycle 

Various multigrid cycles can be used. When developing a multigrid code, 
it is best to use a prescribed cycle so that the results obtained by testing 
different relaxation and interpolation schemes can be easily compared. One 
example of this type is the V-cycle. One V-cycle consists of the following: 
a predefined number of sweeps on each grid In descending order of fineness 
until the coarsest grid Is reached; iterating the coarsest grid problem to 
convergence; and a predefined number of sweeps on each grid in ascending order. 
The goal in multigrid Is to obtain an order of magnitude error reduction per 
cycle. Once the best relaxation and interpolation schemes for the problem have 
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been determined, the adaptive multigrid method can be used. This Is generally 
more efficient than the prescribed cycle method. The adaptive algorithm stays 
on a particular grid until convergence on that grid has slowed to a defined 
"slow" rate, at which time it automatically moves to the next coarser grid. 
Whenever it converges to the set tolerance at a particular grid, it moves back 
up to the next finer grid. 

Determination of Smoothness Between Grids 

The process of determining smoothness between grids Is described using FAS 
and the adaptive multigrid algorithm. Let M denote the finest grid. 

L M (8 M > . F M 

After each relaxation on grid M, the program decides whether to stay on 
the present grid or move to a coarser one. For simplicity sake, this is often 
done by measuring the rate of convergence and determining a cutoff rate, i.e., 
if convergence Is "slow," the program will move to a coarser grid (see 
Fig. 3(a) for a schematic of the general problem). This method Infers that 
slow convergence signifies a smoothed residual function. It assumes that the 
presence of high-frequency terms will show up in a rapidly decreasing global 
error term, where 

global error * 

This Is not a bad assumption If there are no local areas containing 
rapidly changing spatial coefficients. Problems may occur if local regions of 
high-frequency residuals exist within a globally smooth domain. The global 
error term may not be affected by the high-frequency errors and will 
interpolate the solution to the coarser grid, where the residual error of the 
local region will not be well represented. Some other method of determining 
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the smoothness of this local area Is needed in such a case. Additional sweeps 
over the fine grid, or parts of the fine grid, are then needed to smooth the 
error locally. 

Assuming that the residuals have been sufficiently smoothed on the fine 
grid k, the problem moves to grid k-1 and relaxes Eq. (3). The program 
either returns to the fine grid If the solution is converged, or goes to a 
coarser grid if the convergence is slow and the residuals are smooth. When 
the coarsest grid is reached, a converged solution is obtained by continued 
relaxation or by direct solution. When a converged solution is obtained on a 
grid k-1, the solution is updated on grid k using Eq. (4). 

APPLICATION OF MULTIGRID 

The use of implicit time differencing necessitates the solution of an 
NIxMJ set of equations at each time step. In this study, no attempt is made 
to use multigrid across physical time. Multigrid Is used to facilitate the 
iterative solution of the NIxMJ set of equations within each time step. 

The fine grid (M) equations take the form L^{e M } = f m , where L M {6 M } 
represents the differencing scheme of 0 described earlier. The forcing 
function, F m , represents terms from the previous time step which evolve from 
the implicit time differencing scheme. 

F M e*h* 

= At 

The coarse grid representation as derived above are used to implement the 
multigrid procedure. A flow schematic of the procedure used is shown in 
Fig. 3(b). 

Full weighting is used for the fine-to-coarse grid interpolations, taking 
into account all nine fine grid points associated with the coarse grid 
equivalent point, 
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e k 



k 

1-1 J-i 


+ 26 


k 

1-1 J 


+ e 


k 

l-l.j+l 


+ 26 


k 

1 J-l 


+ 46 


k 

1.J 


+ 26 


k 

1,j+1 


+ 6 


k 

1+1 ,j-l 


+ 26 


k 

1 + 1, j 


+ 6 


k 

1+1 ,j+l 




Both linear and third degree polynomial Interpolation from coarse to fine grid 
were tested, and several different relaxation schemes were tried. 

The Switch Function 

The values of the switching function, g(6), are not allowed to change 
during a multigrid solution process. Attempts to let them vary with the 

solution led to major Instabilities. The g values at the new time step, gt, 
are first determined from the previous solution, l.e., 6*“^ v . Using these 
values at time step t, the fine grid equations are relaxed a certain number of 
times, after which the g values are updated to the present ®jpp rox - These 
g values are then used throughout the multigrid solution. 

RESULTS 


Excellent results were obtained over the time steps prior to the start of 
cavitation, both In terms of comparison with a single grid Iterative solution 
and comparison with Brewe's direct solution. These results are summarized In 
Table II. 

Comparisons with single grid Iteration are done on the basis of work units 
(WU) used, where 1 WU is equivalent to 1 relaxation sweep over the finest grid. 
Letting M, i.e., total number of grids, represent the finest grid and k a 
coarser grid, numbered In decreasing order respectively, the equivalent WU 
used by grid k Is 

UUk _ 4<k-M) 


The following results were obtained for the test case having a maximum 
eccentricity of 0.8 and a minimum eccentricity of 0.1 (see Table I). This is 
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one of the more difficult cases to run, since very high pressures and large 
pressure gradients are Induced near the minimum film thickness. The pressure 
gradients and the bubble shape change relatively rapidly with time. 

In all of the trials, a system of three grids was used. The addition of 
a fourth coarsest grid had a negligible effect. For comparison purposes, a 
single grid solution using Gauss-Seldel relaxation on a 96- by 24-point mesh 
was run (96 points circumferentially and 24 points axially). The 96- by 
24-point mesh ensures a grid aspect ratio (AR » Ax/Az) of 1 . It required about 
300 WU per time step. 

Gauss-Seldel (G-S) and Jacobi (J) relaxation schemes with no 
overrelaxatlon were found to be the most effective smoothers for this problem. 
Circumferential line relaxation, l.e., solving simultaneously each line of 
points In the circumferential direction, Is an effective smoother, but Is not 
worth the substantially greater computational time needed to solve for the 
periodic boundary conditions, which Introduce corner terms to the tridiagonal 
matrix. Line relaxation Is used as a local smoother, however, when cavitation 
develops. Both the G-S and J relax points In the direction of the flow, 
i.e., the circumferential direction, sweeping across the axial direction. 
Sweeping across the circumferential direction Is not very effective, nor is a 
combination of axial and circumferential sweeps. A red-black scheme is also 
not very effective. 

The difference between the G-S and the J schemes when used in the 
multigrid process is extremely small. J relaxation uses an average of 0.5 WU 
more than G-S per time step. The reason seems to be that the J multigrid 
uses the same number of fine and medium grid sweeps as does the G-S per 
solution and makes up for its lower efficiency by using a greater number of 
the coarsest grid sweeps, which are very cheap. The advantage of using J 
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when working with a parallel processing computer is that the inner loop is 
vectorizable, since all points in a line are substituted at the end of the line 
instead of point by point. Using the CRAY XMP, one J sweep takes one-tenth 
the CPU time of a G-S sweep. 

A multigrid solution using G-S relaxation along the direction of flow 
and a linear interpolation scheme from the coarse to fine grid was first 
tested. The solution on a 48- by 48-point fine mesh (AR s 3.1) required an 
average of 24 WU per time step. The solution on the 96- by 24-point mesh 
required an average of 14 WU per time step, which is nearly 22 times faster 
than the single grid solution. 

A third degree polynomial interpolation scheme from coarse to fine grid 
was also tested. Using the same 96- by 24-point mesh and G-S, this scheme 
reduces the work per time step to an average of 7.5 WU, half the work used by 
the linear scheme. Also, the third degree polynomial routine takes virtually 
the same amount of CPU time as the linear routine on the CRAY XMP, making it 
highly worthwhile. This scheme used approximately 1/40 of the work used by 
single grid iteration. 

The adaptive multigrid cycle was used to obtain the above results. To 
determine the efficiency, however, a V-cycle was also run. Each cycle reduces 
the error by nearly an order of magnitude. 

The results also compare well with Brewe's direct numerical solution, both 
in accuracy and CPU time. Also, load capacities were compared at various time 
steps. The greatest difference found between the load capacity values is two 
parts in 10 000. Both the direct and the multi grid codes are vectorized to the 
highest efficiency. Both were run on the CRAY XMP for 5000 time steps of 
uncavitated flow. The direct solution took 1086 sec CPU, while the multigrid 
code took 57 sec CPU, about 1/20 the CPU time of the direct solution. 
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The presence of cavltated points In the flow, l.e., the presence of an 
area having g ■ 0 bounded by points having g = 1, requires a more Involved 
approach. On a single grid, the algorithm handles the cavitation as 
efficiently as It would an uncavltated configuration. Problems begin to occur 
when coarser grids are added. 

Initially, the coarse grid cavitation area was determined by injecting 
corresponding fine grid g‘s directly to the coarse grid points. Figure 4 
shows graphs of the residual function at similar states in the solution for 
both an uncavltated and a cavltated configuration using this scheme. These 
graphs were obtained with no extra smoothing around the cavltated area. As can 
be seen, high-frequency local-error terms exist around the cavitated boundary, 
whereas the uncavltated region has already been well smoothed. If the program 
Is allowed to continue from this point, It moves to the coarser grid, where the 
cavitated boundary residuals are not well represented. Depending on how 
unsmooth the boundary residuals are, V-cycle results range from 40-percent 
error reduction per cycle to a slight divergence of error terms per cycle. 

Extra local smoothing helps Immensely, as would be expected. The best 
results are obtained by using a local circumferential line relaxation scheme 
over the cavltated region and boundary points. This scheme Is a very powerful 
smoother and Is also expedient, since, as a local smoother, It reduces to a 
purely tridiagonal matrix of a relatively small number of points. 

The problem still remains of deciding how many local smoothing sweeps is 
"enough," or whether any are necessary at all. If the number of sweeps is set 
such that the most difficult cavitation configurations converge efficiently, 
then configurations having smoother Initial residuals become much less 
efficient. Some type of smoothness Indicator Is necessary. The present 
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routine takes the fine grid residual function, f k , and interpolates it down to 
grid k-1 using the same full weighting as in an actual grid switch. The 
coarse grid values are then interpolated linearly back up to the fine grid k. 
These then represent "smooth" fine grid residuals and are compared with the 
actual residuals. If any of the actual terms fall outside of an envelope 
placed around the smooth term, the problem is deemed unsmooth and local 
smoothing is done. 

The above procedure does much to stabilize the solution process, resulting 
in an average usage of 40 WU per time step solution. For the g-lnjection 
model, however, order of magnitude error reduction per cycle is not usually 
obtained, and certain cavitation configurations do occur which are very 
difficult or Impossible to solve. 

As mentioned earlier, the coarse-grid operator need not be the same 

as the fine-grid operator lA Because of this fact, some latitude in handling 
the values of the g coefficient In the coarse grid equations is permissible. 
This led to a coarse-grid determination scheme for the g values that not only 
circumvented occurrences of Instability due to Injection but had a major 
beneficial effect on the. solution. Recall, the g values on the finest grid 
are determined by the value of 6 at each nodal point, i.e., g has a value of 
1 at full film points and a value of 0 at cavltated points. It was found that 
stability of the solution across all possible cavitation configurations can be 
obtained by defining a parameter FG as: 


FG 


' k k k k „k k k 

= [ g i-l,j-l + g i-l,j + 9 1-l,j+l + 9 i,j-l + 9 1,j + g i , j+1 + g 1+l,j-l 


+ 9 i+i,j + 9 


k 1 

1+1 ,3+1 J 


If 

FG = 0; 

g k -' - 0 

If 

FG * 0; 

- 1 
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In other words, a fine grid point must have a g k value of 0 and must be 
surrounded < al 1 ' eight points) by points having g k = 0 in order for the 
corresponding coarse grid point to be set to g k_1 = 0. Other schemes were 
found to work but were not as efficient. This scheme resulted in an average of 
20 WU per time step, the number of WU ranging from 11 to 35 WU. The solution 
process remains stable throughout bubble formation and bubble collapse. The 
time steps which required the most work units occurred at the very beginning of 
bubble formation, when there were very few cavltated points, and during bubble 
collapse. While the bubble is collapsing, it Is also experiencing its greatest 
amount of movement downstream, so It might be this movement rather than the 
process of collapse which requires more solution time. 

A V-cycle analysis shows that better than an order of magnitude error 
reduction per cycle is obtained, though at the cost of extra smoothing on the 
finer grids. 

The results for cavltated flow also compare well with Brewe's direct 
solution. The same bubble shape, motion, and duration are obtained from both 
programs. Figure 5 contains computed pressure distributions at various time 
steps. The cavltated area Is Indicated by the flat area of zero pressure 
gradient. Comparison of load capacity terms shows a maximum difference of five 
parts in 10 000 In the cavitated region. When run on the CRAY XMP, the direct 
solution of 5000 time steps of cavltated flow again takes 1086 sec CPU. The 
multigrid solution of 5000 time steps of cavltated flow takes 150 sec CPU, 
approximately one-eighth the CPU taken by the direct solution. Even though the 
cavitated configurations take more CPU than do the uncavitated configurations, 
the multigrid solution still represents a very significant savings over the 
direct method. 
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In his paper on the direct ADI solution of the Elrod algorithm, Brewe (2) 
states that his solution uses two to three times the computational time used 
by an iterative solution of a nonconservative film model using pseudo-Giimbel 
boundary conditions. This nonconservative film model is only suited to steady- 
state conditions, but is often used by industry. Thus, the multigrid solution 
of the Elrod algorithm requires about one-tenth to one-third the computational 
time (for uncavitated and cavitated flow respectively) of the nonconservative 
film model solution, while still retaining the more realistic representation of 
the flow. 

CONCLUSION 

A multigrid iterative technique is used in the solution of the Eirod 
algorithm for the case of a dynamically loaded journal bearing undergoing 
cavitation. This solution is compared both to a single grid iterative solution 
in terms of work used, and to a direct ADI solution in terms of computer time 
required. Excellent results are obtained both prior to and during cavitation, 
although the presence of cavitation does introduce difficulties in the solution 
process. 

The best results are obtained using the following: a grid aspect ratio of 

1; full weighting interpolation from the fine grid to the coarse grid; third 
degree polynomial interpolation from the coarse grid to the fine grid; either 
Gauss-Seidel or Jacobi relaxation with no overrelaxation. Implementing these 
techniques, the solution at time steps prior to cavitation uses 1/40 the amount 
of work used by a single grid iterative Gauss-Seidel solution and 1/20 the 
computer time used by the direct ADI solution. During cavitation, the 
multigrid solution uses 1/15 the amount of work used by a single grid G-S 
solution and one-eighth the computation time used by the direct ADI solution. 
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Based on the results stated In this paper, it is evident that the solution 
of the Elrod algorithm using multigrid techniques provides an extremely viable 
method to industry for the solution of journal bearing problems. 
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TABLE I. - OPERATING CONDITIONS 


AR, m . . . . 

. . 5.0xl0~ 4 

R, m . . . . 

. . . 0.0425 

L/D 

.... 1.0 

e 

. 0.1 to 0.8 

u s , rad/s . . 

. . . -19.5 

con , rad/s . . 

... -92.7 

13, N/m 2 . . . 

. . 1 .72x1 0^ 

p, N-s/m 2 . . 

. . . 0.066 

p a , N/m 2 . . 

. 1.0133x10 s 

p c , N/m 2 . . 

.... 0.0 
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TABLE II. - RESULTS 


(a) Prior to cavitation, work units (WU a ) used per time step 
solution. 


Number 

of 

grids 

Aspect 

ratio, 

AR 

Type of 
relaxation 

Type of 
interpolation 
k to k+1 

Average 

WU 

1 

l 



300 

24 

14 

7.5 

8.0 

3 

3.1 

1 

1 

1 

VJ J 

J 

Linear 

Linear 

Third-degree 

polynomial 

Third-degree 

polynomial 


a l WU = the equivalent to one relaxation sweep over a 96- by 
24- (or a 48- by 48-) point grid. 


(b) Prior to cavitation, CPU time used for 
solution of 500 time steps on Cray XMP. 


Type of solution 

CPU 

time, 

sec 

Direct - ADI (96- by 24-point mesh) 

1086 

Multigrid - three grids (Jacobi 

57 

relaxation, third-degree poly- 


nomi a 1 i nterpol at i on from coarse 


to fine, 96- by 24-point mesh) 



(c) During cavitation 


Number 

of 

grids 

Type of 
solution 

g-model b 

Average 
WU per 
time step 
solution 

Cray CPU 
time for 
5000 time 
steps, 
sec 

Stabi 1 i ty of 
solution 
process 



G-S 

g = fcn(0) 

300 

— 

Stable 


3 

G-S or J; 

Injection 


— 

Unstable 



no L c 








G-S or J 

Injection 

40 

— 

Unstable 



and L 








G-S or J 

gi 

40 to 50 

— 

Stable 



and L 








G-S or J 

gsur 

20 

150 





and L 






1 

Direct (ADI) 

g = fcn(6) 


1086 



1 

Nonconformal 

Stationary 


N500 





film modi- 

cavita- 







fication 

tion area 






b G-mode1s : 

Injection: q M = fcn(0) = 1 or 0; g^< M = corresponding g^ + '. 

gl : g M = fcn(6); all g k<M = 1. 

gsur: g M = fcn(0); g^<M = 0 only if corresponding g k +l = 0 and is 

k+1 

surrounded by g =0 points. 

C L = local circumferential line relaxation around cavitated area. 







CONTROL VOLUME CV AT POINT X 
FIGURE 2. - CONTROL VOLUME APPROACH. 





(b) MULTIGRID-ELROD ALGORITHM SCHEMATIC (BOUNDARY CONDITIONS 
NOT INCLUDED). 


FIGURE 3. - MULTIGRID ALGORITHM SCHEMATICS WITH FULL APPROXIMATION STORAGE. 
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FIGURE 5. - PRESSURE DISTRIBUTION AT VARIOUS TIME STEPS. 
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